Non— equilibrium Anisotropic Phases, Nucleation and Critical Behavior in a Driven 

Lennard— Jones Fluid 
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We describe short-time kinetic and steady-state properties of the non-equilibrium phases, namely, 
solid, liquid and gas anisotropic phases in a driven Lennard-Jones fluid. This is a computationally- 
convenient two-dimensional model which exhibits a net current and striped structures at low temper- 
ature, thus resembling many situations in nature. We here focus on both critical behavior and details 
of the nucleation process. In spite of the anisotropy of the late-time "spinodal decomposition" pro- 
cess, earlier nucleation seems to proceed by Smoluchowski coagulation and Ostwald ripening, which 
are known to account for nucleation in equilibrium, isotropic lattice systems and actual fluids. On 
the other hand, a detailed analysis of the system critical behavior rises some intriguing questions on 
the role of symmetries; this concerns the computer and field-theoretical modeling of non-equilibrium 
fluids. 



I. INTRODUCTION 

Steady states in non-equilibrium many-particle sys- 
tems typically involve a constant flux of matter, charge, 
or some other quantity and, consequently, stripes 
or other spatial anisotropics show up at appropriate 
scales^i^i^. This occurs during segregation in driven 
sheared systems^i^, flowing fluidsS', shaken granular 
matteriSiiS, and non-equilibrium liquid-liquid binary 
mixturesii, and it has been reproduced in computer sim- 
ulations of driven colloida]~ and fluid^'^^ systems, for 
instance. Further examples are the anisotropics observed 
in both high-temperature superconductorsi^si^ and elec- 
tron gasesiSiii. The ripples shaped by the wind in sand 
desertsiSiiS and the lanes and trails formed by living or- 
ganisms and vehicle trafRo22iSi also share some of the 
essential physics. 

Lacking theory for the "thermodynamic" instabilities 
causing the observed striped structures, one tries to link 
them to the microscopic dynamics of suitable model sys- 
tems. For two decades, the driven lattice gas (DLG)S^, 
namely, a computationally-convenient model system 
in which particles diffuse under an external driving 
"field" , has been a theoretical prototype of anisotropic 
behaviori*22iS4. This model was recently shown to be 
unrealistic in some essential sense, however'^^. Particle 
moves in the DLG are along the principal lattice direc- 
tions, and any site can hold one particle at most, so that 
a particle impedes the one behind to jump freely along 
the direction which is favored to model the action of the 
field. Consequently, the lattice geometry acts more effi- 
ciently in the DLG as an ordering agent than the field 
itself, which occurs rarely — never so dramatically — in 
actual cooperative transport. In fact, actual situations 
may in principle be more closely modeled by means of 
continuum models, and this peculiarity of the DLG im- 
plies that it lacks a natural off-lattice cxtcnsion^^. 

Here we present, and analyze numerically a novel non- 
equilibrium off-lattice, Lennard-Jones (LJ) system which 



is a candidate to portray some of the anisotropic behav- 
ior in nature. The model, which involves a driving field 
of intensity E, reduces to the celebrated (equilibrium) 
LJ flm(&i^ as E -^ 0. For any E > 0, however, it 
exhibits currents and anisotropic phases as in many ob- 
servations out of equilibrium. In particular, as the DLG, 
our model in two dimensions shows striped steady states 
below a critical point. We also observe critical behavior 
consistent with the equilibrium universality class. This 
is rather unexpected in view of the criticality reported 
both for the DLG and in a related experimenlj^. On the 
other hand, concerning the early-time relaxation before 
well-defined stripes form by spinodal decomposition, we 
first observe — as in previous studies of relaxation to- 
wards equilibrium — effective diffusion of small droplets, 
which is followed by monatomic diffusion probably com- 
peting with more complex processes. It is very likely that 
our observations here concerning nucleation, coexistence, 
criticality, and phases morphology hold also in a number 
of actual systems. 

The paper is organized as follows. In section II we de- 
fine the model, and section III is devoted to the main 
results as follows. § III. A describes the early-time segre- 
gation process as monitored by the excess energy, which 
measures the droplets surface. § III.B describes some 
structural properties of the steady state, namely, the ra- 
dial and azimuthal distribution functions, and the de- 
gree of anisotropy. § III.C, which depicts some trans- 
port properties, is devoted to an accurate estimate of the 
liquid-vapor coexistence curve and the associated critical 
indexes. Section IV contains a brief conclusion. 



II. THE MODEL 

Consider N particles of equal mass (set henceforth to 
unity) in a d— dimensional box, L'^, with periodic bound- 
ary conditions. Interactions are via the truncated and 
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FIG. 1: Schematic representation of the region (grey) which 
is accesible to a given particle as a consequence of a trial move 
for _B = (left) and E = oo (right), assuming the "infinite" 
field points x, horizontally. 



FIG. 2: Typical steady-state configurations for E = (top 
row) and E ^ oo (bottom row) at T = 0.10, 0.15, 0.30 and 
0.35, respectively, from left to right. This is for A*' = 1000 
and p — 0.30. 



shifted pair potentialSi: 



where 



(l)Lj{r) - (j>Lj{rc) if r < Tc 
if r ^ Tc, 



0L,7(r) = 46 [(a/r)i2 - {a/rf 



(1) 
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is the LJ potential, r stands for the interparticle distance, 
and Tc is a cut-off that we set Vc = 2.5a. The parameters 
a and e are, respectively, the characteristic length and 
energy — that we use in the following to reduce units as 
usual. 

Time evolution is by microscopic dynamics according 
to the transition probability per unit time (rate): 

J^) [^ ^ ^,) = ;^(^) X min {l, e-^*'/^} . (3) 

Here, 
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E is the intensity of a uniform external field along a prin- 
cipal lattice direction, say x, 1] = {ri, . . . , r/v} stands for 
any configuration of energy 
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where r^ is the position of particle i that can be anywhere 
in the d— torus, rji equals rj except for the displacement of 
a single particle by Si — r[~ ri, and A$i = $(r/i) — ^{i]) 
is the cost of such displacement. 

It is to be remarked that x'"^\ ^s defined in 10}, con- 
tains a drive bias (see Fig.QJ such that the rate Q lacks 
invariance under the elementary transitions rj ^ rji. Con- 
sequently, unlike in equilibrium, there is no detailed bal- 
ance for toroidal boundary conditions if _E > 0. 

We report here on the results from a series of 
Monte Carlo (MC) simulations using a neighbor-list 
algorithm.^''. Simulations concern fixed values of N, with 
N < 10'*, particle density p = N/L'^ within the range 



p G [0.2,0.6], and temperature T G [lO-^ 10^] . Fol- 
lowing the fact that most studies of striped structures, 
e.g., many of the ones mentioned in the first paragraph 
of section ^ concern two dimensions — in particular, the 
DLG critical behavior is only known with some confi- 
dence for d = 34j28j29 — .^g restricted ourselves to a two 
dimensional torus. The maximum particle displacement 
is (5max = 0.5 in our simulations. We report below on 
steady-state averages over 10^ configurations, and ki- 
netic or time averages over 40 or more independent runs. 

The distribution of displacements Si is uniform, except 
that the new particle position f^ is (most often in our 
simulations) sampled only from within the half-forward 
semi-circle of radius ^max centered at fi , as illustrated in 
the right graph of Fig. ^ This is because the infinite- 
field limit, E -^ cxD, turns out to be most relevant, and 
this means, in practice, that any displacement contrary 
to the field is forbidden. This choice eliminates from the 
analysis one parameter and, more importantly, it hap- 
pens to match a physically relevant case. As a matter 
of fact, simulations reveal that any external field E > 
induces a flux of particles along x — which crosses the 
system with toroidal boundary conditions — that mono- 
tonically increases with E, and eventually saturates to 
a maximum. This is a realistic stationary condition in 
which the thermal bath absorbs the excess of energy dis- 
sipated by the drive. 



III. MAIN RESULTS 

Fig.|21illustrates late-time configurations, i.e., the ones 
that typically characterize the steady state, as the tem- 
perature T is varied. These graphs already suggest that 
the system undergoes an order-disorder phase transition 
at some temperature Te- This happens to be of second 
order for any £' > 0, as in the equilibrium case E — 0. 
We also observe that Te decreases monotonically with in- 
creasing E, and that it reaches a well-defined minimum. 

Too, SiS E —>■ GO. 

Fig. 121 also shows that, at low enough temperature, an 




FIG. 3: Typical configurations for E — (top row) and 
E ^ oo (bottom row) as time proceeds during relaxation 
from a disordered state (as for T —> oo) at i = 0. The graphs 
are, respectively from left to right, for times t — 10^, 10*, 
10^ and 10'' MC steps. This is for TV = 7500, p = 0.35, and 
T = 0.275, below the corresponding transition temperature. 



anisotropic interface forms between the condensed phase 
and its vapor; this extends along x throughout the system 
at intermediate densities. 




FIG. 4: Time evolution of the enthalpy per particle for A'^ = 
7500, p = 0.35 and, from top to bottom, T = 0.200, 0.225, 
0.250 and 0.275. Straight lines are a guide to the eye; the 
slope of each line is indicated. The inset shows the detail 
at early times. (For clarity of presentation, the main graph 
includes a rescaling of the time corresponding to the data for 
T = 0.250, 0.225 and 0.200 by factors 2, 3 and 3, respectively.) 



A. Phase segregation kinetics 

Skipping microscopic details, the kinetics of phase seg- 
regation at late times looks qualitatively similar to the 
one in other non-equilibrium cases, including driven lat- 
tice systemsi^ and both molecular-dynamic'^'' and Cahn- 
Hilliard^ representations of sheared fluids, while it essen- 
tially differs from the one in the corresponding equilib- 
rium system. This is illustrated in Fig. |S1 One observes, 
in particular, condensation of many stripes — as in the 
graph for t = 10^ in Fig. — into a single one — as in 
the first three graphs at the bottom row in Fig. |2| This 
process corresponds to an anisotropic version of the so- 
called spinodal decomposition^^, which is mainly charac- 
terized by a tendency towards minimizing the interface 
surface as well as by the existence of a unique relevant 
length, e.g., the stripe widtbi^. A detailed analysis of 
this late regime, which has already been studied for both 
equilibriumSi^ and non-equilibrium cases, including the 
j^LQi3j35^ will be the subject of a separate report. 

Detailed descriptions of early non-equilibrium nucle- 
ation are rare as compared to studies of the segrega- 
tion process near completion. Following an instantaneous 
quench from a disordered state into T < Too {p) , one ob- 
serves in our case that small clusters form, and then some 
grow at the expenses of the smaller ones but rather inde- 
pendently of the growth of other clusters of comparable 
size. This corresponds to times t < 10^ in Fig. O i.e., 
before many well-defined stripes form. We monitored in 
this regime the excess energy or enthalpy, H (t) , mea- 
sured as the difference between the averaged internal en- 
ergy at time t > and its stationary value. This reflects 
more accurately the growth of the condensed droplets 
than its size or radius, which are difhcult to be estimated 
during the early stagea^Si2i. Furthermore, H (t) may be 
determined in microcalorimetric experiments^. 



The time development of the enthalpy density h (t) — 
H (t) /N is depicted in Fig. ^ This reveals some well- 
defined regimes at early times. 

The first regime, (a) in the inset of Fig. ^ follows a 
power law t~^ with 9 « 0.165 — which corresponds to the 
line shown in the graph — independently of the temper- 
ature investigated. This is the behavior predicted by the 
Smoluchowski coagulation or effective cluster diffusionS^. 
The same behavior was observed in computer simulations 
for £^ = and also reported to hold in actual experi- 
ments on binary mixtures^S*^. This suggests the early 
dominance of a rather stochastic mechanism, in which 
the small clusters rapidly nucleate, which is practically 
independent of the field, i.e., it is not affected in practice 
by the drive. The indication of some temperature depen- 
dence in equilibriurr^, which is not evident here, might 
correspond to the distinction between deep and shallow 
quenches made in Ref,^ that we have not investigated 
out of equilibrium. 

At latter times, there is a second regime, (b) in Fig. 01 
in which the anisotropic clusters merge into filaments 
and, finally, stripes. We observe in this regime that 
9 varies between 0.3 and 0.6 with increasing T. Ost- 
wald ripening^, consisting of monomers diffusion, pre- 
dicts 9 = 1/3. It is likely that regime (b) describes 
a crossover from a situation which is dominated by 
monomers at low enough temperature to the emergence 
of other mechanisms^^i^ which might be competing as 
T is increased. 

Finally, one observes a regime, (c) in Fig. ^ which 
corresponds to the beginning of spinodal decomposition. 
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FIG. 5: Details of the structure in the low—T, solid phase as 
obtained by zooming in into configurations such as the ones 
in Fig. H This is for T = 0.05, 0.12 and 0.25, from left to 
right, respectively. 



B. Structure of the steady state 

For any E > i)^ the anisotropic condensate changes 
from a solid-hke hexagonal packing of particles at low 
temperature (e.g., T = 0.10 in Fig. EJ, to a poly- 
crystalline or perhaps glass-like structure with domains 
which show a varied morphology at (e.g.) T — 0.12. The 
latter phase further transforms, with increasing temper- 
ature, into a fluid-like structure at (e.g.) T = 0.30 and, 
finally, into a disordered, gaseous state. 

More specifically, the typical situation we observe at 
low temperature is illustrated in Fig. O At sufhciently 
low temperature, T = 0.05 in the example, the whole 
condensed phase orders according to a perfect hexagon 
with one of its main directions along the field direction 
X. This is observed in approximately 90% of the config- 
urations that we generated at T = 0.05, while all the 
hexagon axis are slanted with respect to x in the other 
10% cases. As the system is heated up, the stripe looks 
still solid at T = 0.12 but, as illustrated by the second 
graph in Fig.[Sl one observes in this case several coexist- 
ing hexagonal domains with different orientations. The 
separation between domains is by line defects and/or va- 
cancies. Interesting enough, as it will be shown later on, 
both the system energy and the particle current are prac- 
tically independent of temperature up to, say T ~ 0.12. 
The hexagonal ordering finally disappears in the third 
graph of Fig. [SJ which is for T = 0.25; this case corre- 
sponds to a fluid phase according to the criterion below. 

A close look to the structure is provided by the radial 
distribution (RD), 
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i.e., the probability of finding a pair of particles a distance 
r apart, relative to the case of a random spatial distri- 
bution at same density. This is shown in the lower inset 
of Fig. El At fixed T, the driven fluid is less structured 
that its equilibrium counterpart, suggesting that the field 
favors disorder. This is already evident in Fig. |21 and it 
also follows from the fact that the critical temperature 
decreases with increasing E. 

The essential anisotropy of the problem is revealed by 




FIG. 6: Data from simulations for TV = 7000 and p = 0.35. 
The main graph shows the degree of anisotropy, as defined in 
the main text, versus temperature. The vertical dotted line 
denotes the transition temperature. The lower (upper) inset 
shows the radial (azimuthal) distribution at T = 0.20, full 
line, and T — 0.30, dashed line. 



the azimuthal distribution (AD) defined 



\i<J 



(7) 



where 6*^ G [0, 2tt) is the angle between the line con- 
necting particles i and j and the field direction x. Ex- 
cept at equilibrium, where this is uniform, the AD is 
7r/2— periodic with maxima at kir and minima at fc7r/2, 
where k is an integer. The AD is depicted in the upper 
inset of Fig. ^ 

We also monitored the degree of anisotropy, defined as 
the distance 
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which measures the deviation from the equilibrium, 
isotropic case, for which a{9) ~ 1, independent of 9. 
The function ||SJ), which is depicted in the main graph 
of Fig. |S1 reveals the existence of anisotropy even above 
the transition temperature. This shows the persistence of 
non-trivial two-point correlations at high temperatures 
which has been demonstrated for other non-equilibrium 
models^. 



C. Coexistence curve 

The transition points may be estimated from the tem- 
perature dependence of the mean potential energy per 
particle. 



/ = 7V-1 ($(,,)) , 



(9) 



and from the net current j, defined as the mean displace- 
ment per MC step per particle. Fig.[7|shows well-defined 




FIG. 7: Temperature dependence of the mean energy 
(squares; the scale is on the right axis) and normalized net 
current (circles; scale at the left) for N — 5000 and p = 0.30 
under "infinite" field. The inset shows the T— dependence of 
the current over a wider range. 



changes of slope in both magnitudes when the phase 
transforms from sohd to hquid (T « 0.12) and then to 
disorder {T « 0.30). The persistence of correlations is 
again revealed by the fact that the current is nonzero for 
any, even low T, though it is small, and roughly indepen- 
dent of T, in the solid phase. The energy Q behaves lin- 
early with temperature for T E (0.12,0.3) , as expected 
for a fluid phase. The maximum value of the current, 
jmax = ^^max/Sn, is Only rcachcd for T —>■ cxd. The way 
this limit is approached is illustrated in the inset of Fig. [3 
where the grow is shown to be slower than exponential. 

A main issue concerning the steady state is the liquid- 
vapor coexistence curve and the associated critical behav- 
ior. The (non-equilibrium) coexistence curve may be de- 
termined from the density profile transverse to the field. 
This is illustrated in Fig. |H1 

At high enough temperature — in fact, already at 
T = 0.35 in this case for which the transition temper- 
ature is slightly above 0.3 — the local density is roughly 
constant around the mean system density, p = 0.35 in 
Fig. |S1 As T is lowered, the profile accurately describes 
the existence of a single stripe of condensed phase of den- 
sity p+ which coexists with its vapor of density p_ . The 
interface becomes thinner and smother, and p+ increases 
while p- decreases, as T is decreased. 

As in equilibrium, one may use p-^- — p- as an order 
parameter. The result of plotting p+ and p- at each tem- 
perature results in the non-symmetric liquid-vapor coex- 
istence curve shown in Fig. |51 The same result follows 
from the current, which in fact varies strongly correlated 
with the local density. Notice that the accuracy of our 
estimate of p± is favored by the existence of a linear in- 
terface. This is remarkable because we can therefore get 
closer to the critical point than in equilibrium. Further- 
more, we found that the rectilinear diameter law, 
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FIG. 8: Density profiles transverse to the field for A'' = 7000, 
p = 0.35, and diff'erent temperature, as indicated. The coex- 
isting densities, p±, are indicated. 



expansio: 
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and the scaling law (the first term of a Wegner-type 



can be used here to estimate the critical parameters with 
higher accuracy than in the equilibrium case^. The sim- 
ulation data in Fig. O then yields the values in Table I, 
which are confirmed by the familiar log-log plots. Com- 
pared to the equilibrium critical temperature reported 
by Smit and Frenkel2&, one has that Tq/T^o ~ 1.46, i.e., 
the change is opposite to the one for the DLGA. This 
confirms the observation above that the field acts in the 
non-equilibrium LJ system favoring disorder. 

The fact that the order-parameter critical exponent 
is relatively small may already be guessed by noticing 
the extremely flat coexistence curve in Fig. |51 This is 
similar to the corresponding curve for the equilibrium 
two-dimensional LJ fluidsSSiiSiH, and it is fully consis- 
tent with the equilibrium Onsager value, fi = 1/8. We 
therefore believe that our model belongs to the Ising uni- 
versality class. In any case, one may discard with confi- 
dence both the DLG value (5 k, 1/3 as well as the mean 
field value (3 — 1/2 which was reported for fluids under 
shear— — both cases would produce a hump visible to the 
naked eye in a plot such as the one in Fig. |51 One may 
argue that this result is counterintuitive, as our model 
apparently has the short-range interactions and symme- 
tries that are believed to characterize the DLG. 



IV. CONCLUSION 

In summary, the present (non-equilibrium) two- 
dimensional Lennard-Jones system, in which particles 
are subject to a constant driving field, has two main 
general features. On one hand, this case is more con- 
venient for computational purposses, than others such 
as, for instance, standard molecular-dynamics realiza- 
tions of driven fluid systems. On the other, it seems to 
contain the necessary essential physics to be useful as a 
prototypical model for anisotropic behavior in nature. 
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FIG. 9: Coexistence curve (squares) for the LJ non- 
equilibrium model obtained from the density profiles in Fig. |S| 
The fluid phase and the coexistence region are indicated. The 
triangles are the arithmetic mean points, which serve to com- 
pute the critical parameters. The large circle at the top of 
the curve locates the critical point, and the solid line is a fit 
using the Wegner expansion and the rectilinear diameter law 
with the critical parameters given in Table I. 



This model reduces to the familiar LJ case for zero 
field. Otherwise, it exhibits some arresting behavior, in- 
cluding currents and striped patterns. We have identified 
two processes which seem to dominate early nucleation 
before anisotropic spinodal decomposition sets in. In- 
teresting enough, they seem to be identical to the ones 



characterizing a similar situation in equilibrium. 

We have also concluded that the model critical behav- 
ior is consistent with the Ising one for d — 2 but not 
with the critical behavior of the driven lattice gas. This 
is puzzling. For instance, using the language of statis- 
tical field theory, symmetries seem to bring our system 
closer to the non-equilibrium lattice model than to the 
corresponding equilibrium case. The additional freedom 
of the present, off-lattice system, which in particular im- 
plies that the particle-hole symmetry is violated — which 
induces the coexistence-curve asymmetry in Fig.|51in ac- 
cordance with actual systems — are likely to matter more 
than suggested by some naive intuition. 

Further study of the present non-equilibrium LJ sys- 
tem and its possible variations is suggested. A principal 
issue to be investigated is the apparent fact that the full 
non-equilibrium situations of interest can be described 
by some rather straightforward extension of equilibrium 
theory. We here report on some indications of this con- 
cerning early nucleation and properties of the coexistence 
curve. No doubt it would be interesting to compare 
more systematically the behavior of models against the 
varied phenomenology which was already reported for 
anisotropic fiuids. This should also help a better under- 
standing of non-equilibrium critical phenomena. 

We acknowledge very useful discussions with M. A. 
Munoz and F. de los Santos, and financial support from 
MEyC and FEDER (project FIS2005-00791). 
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